package gbench.sandbox.matlib;

import org.junit.jupiter.api.Test;

import java.util.concurrent.atomic.AtomicInteger;

import gbench.common.matlib.MatlibCanvas;
import gbench.common.matlib.MatlibCanvas.NVec;
import gbench.common.matlib.MatlibCanvas.Range;
import gbench.common.matlib.MatlibCanvas.XVec;
import gbench.common.matlib.matrix.Complex;
import gbench.common.matlib.matrix.IComplex;
import gbench.common.matlib.matrix.ComplexMatrix;
import gbench.common.matlib.matrix.IComplex.Relation;

import static java.lang.Math.*;
import static gbench.common.matlib.MatlibCanvas.println;
import static gbench.common.matlib.MatlibCanvas.Range.rg;
import static gbench.common.matlib.matrix.Complex.C;
import static gbench.common.matlib.matrix.Complex.C1;
import static gbench.common.matlib.matrix.Complex.C2;
import static gbench.common.matlib.matrix.Complex.EXP;
import static gbench.common.matlib.matrix.Complex.IMGONE;
import static gbench.common.matlib.matrix.Complex.O;
import static gbench.common.matlib.matrix.Complex.ONE;
import static gbench.common.matlib.matrix.Complex.UNIT;
import static gbench.common.matlib.matrix.ComplexMatrix.KJX;
import static gbench.common.matlib.matrix.ComplexMatrix.V;


/**
 * 复数计算
 * @author gbench
 *
 */
public class JunitComplex {
    
    /**
     * 傅里叶变换
     */
    @Test
    public void fourier() {
        final var vx = NVec.of(1000,i->i);
        final var n = vx.size();
        final var fX = Complex.fourierX(n);
        final var a = fX.mmult(vx); // 傅里叶变换
        final var fY = Complex.fourierX(n, true);
        println("fY\n", fY.mmult(a));
    }
    
    /**
     * 傅里叶变换2
     * omega.mul(k*j).exp() <==> omega.exp().pow(k*j)
     */
    @Test
    public void fourier2(){
        final var n = 10;
        final var omega = C2(2*PI/n);
        final var x = NVec.of(n,i->i);
        println("x",x);
        final var ft = XVec.OF(n, j-> XVec.OF(n, k->omega.mul(k*j).exp().mul(x.get(k)).scale(1/sqrt(n))).reduce2(IComplex::add));
        println("傅里叶变换 ft",ft);
        final var ift = XVec.OF(n, j-> XVec.OF(n, k->omega.exp().pow(-k*j).mul(ft.get(k)).scale(1/sqrt(n))).reduce2(IComplex::add));
        println("逆傅里叶变换 ift",ift);
    
    }

    /**
     * 基本数据运算(代数运算)
     */
    @Test
    public void algebra() {
        final var mx = ComplexMatrix.of(2,3, 
            1,2,3,
            4,5,6);
        println(mx);
        println("add----------------\n",mx.add(mx));
        println("sub----------------\n",mx.sub(mx));
        println("mul----------------\n",mx.mul(mx));
        println("div----------------\n",mx.div(mx));
        println("div mul----------------\n",mx.div(mx).mul(mx));
        println("mul----------------\n",mx.mul(C(1,4)));
        println("add----------------\n",mx.add(C(1,4)));
        println("mmult----------------\n",V(0,1,2,3,4,5,6).mmult(V(0,1,2,3,4,5,6).transpose()));
        println("add----------------\n",C(1,2).add(C(3,4)));
        println("rotate----------------\n",C1(2).rotate(PI/2));
        println("rotate----------------\n",C2(2).rotate(-PI/2));
        println("scale,rotate----------------\n",C2(2).scale(100).rotate(-PI/2));
        println("affine",O.add(C(1).affine(1000,PI/2)));
    }
    
    /**
     * 数据zip操作
     */
    @Test
    public void zip() {
        final var xx = NVec.of(100,i->i+1);
        xx.zipS(xx).map(Complex::new).forEach(e->{
            println(e.affine(100, PI/9));
        });
        
        Complex.C(1).nthRoot(10).forEach(System.out::println);
    }
    
    /**
     * 求方程根
     */
    @Test
    public void root() {
        final var ai = new AtomicInteger();
        Complex.C(8).nthRoot(3).forEach(e->{
            println(ai.getAndIncrement(),e,e.getArgument(),e.theta());
        });
    }
    
    /**
     * 矩阵数据提取
     */
    @Test
    public void matrix() {
        final var mx = ComplexMatrix.of(3,3, 
            1,2,3,
            4,5,6,
            7,8,9);
        println("mx\n",mx);
        println("mx.cgets(1)\n",mx.cgets(1));
        println("mx.rgets(1)\n",mx.rgets(1));
        println("mx.rgets2(1)\n",mx.rgets2(1));
        println("mx.cgets2(1)\n",mx.cgets2(1));
        println("mx.minorX(1, 1)\n",mx.minorX(1, 1));
        println("mx.minorX(1, 1).minorX(0, 0)\n",mx.minorX(1, 1).minorX(0, 0));
        println("mx.minorX(1, 1).minorX(0, 0).minorX(0, 0)\n",mx.minorX(1, 1).minorX(0, 0).minorX(0, 0));
        println("\n----------余子式minors--------------");
        Range.cph(3,3).forEach(aa->{
            println(NVec.of(aa)+"\n",mx.minorX(aa[0], aa[1]));
        });
        println("行列式求值 det\n", ComplexMatrix.of(2, 2, 1, 2, 3, 4).det());
    }
    
    /**
     * 复数绘图
     * @throws InterruptedException
     */
    @Test
    public void plot() throws InterruptedException {
        final var size = 1000;
        final var matlib = new MatlibCanvas(size);
        matlib.plot(size, size / 2, sess -> {
            final var n = 1000;
            final var omega = Math.PI * 2 / n; // 本原单位根的幅角主值
            XVec.OF(n, j -> C2(omega * j).exp().scale(100).conjugate())
                    .map(c -> NVec.of(c.real(), c.imaginary(), 255, 0, 0).affine_offset(size / 2, size / 4))
                    .forEach(sess::drawPoint);

        });// plot

        Thread.sleep(5_000);
    }
    
    /**
     * 简单傅里叶矩阵
     */
    @Test
    public void bar() {
        final var v = V(1,0,-1,0); // 数据向量
        final var n = v.nrows(); // 向量的维度
        final var w = EXP(2*PI/n); // 本原n词方根
        final var Fn = KJX(4, 4, (j,k)->w.pow(j*k)).div(Math.sqrt(n));
        println("傅里叶变换\n",Fn.mmult(v));
        
        final var alpha = V(Range.rg(n).map(Complex::C1)); // 基本空间向量
        final var zeta = alpha.mmult(alpha.transpose()); // 指数矩阵
        println("\npn\n",zeta);
        println("\n傅里叶系数矩阵:jks\n",KJX(n,n,(j,k)->w).pow(zeta));
        println("\n傅里叶系数矩阵:Pn\n",zeta.pow2(w));
    }
    
    /**
     * 接口测试
     */
    @Test
    public void f1() {
        final var a = V(1, 2, 3, 4, 5).dot(V(2, 3, 4, 5, 6));
        println("\n---------1---------\n",a);
        println("\n---------2---------\n",a.oneXone().mul(V(1, 2, 3, 4)).minorX(0, 0).isEmpty());
        println("\n---------3---------\n",a.oneXone().zipS(a.oneXone(), (x, p) -> x));
        println("\n---------4---------\n",a.repeatV(100).T().fmap(d -> d.add(1)));
        final var one = Complex.ONE;
        println("\n---------5---------\n",one.right(5));
        one.repeatV(100);
        println("\n---------6---------\n",one.repeatV(100).add(1).transpose());
        println("\n---------7---------\n",ONE.fx(IMGONE::add).equals(UNIT));
        println("\n---------8---------\n",ONE.swap());

        println("\n---------9---------\n",ONE.compareTo(ONE.add(1)));
        println("\n---------10---------\n",ONE.mul(2).rotate(PI / 2));
        println("\n1---------1---------\n",ONE.relative(ONE.mul(2).rotate(PI / 2)));
        println("\n---------12---------\n",Relation.TP.compareTo(Relation.BT));
        final var d = V(1, 2, 3, 4, 5, 6).fmap(IComplex::sqrt);
        println("\n---------13---------\n",d.dot(d));
        println("\n---------14---------\n",d.mmult(d.tp()));
    }

    /**
     * 接口测试
     */
    @Test
    public void f2() {
        final var a = V(1,2,3,4,5,6);
        final var b = a.fmap(IComplex::sqrt);
        final var c = a.fmap(C(2)::pow);
        final var d = V(Range.rg(-PI,PI).limit(1000)).fmap(IComplex::cos);
        println("\n---------15---------\n",b,c,d);
        println("\n---------16---------\n",V(rg(1,10)).fmap(IComplex::reciprocal));
        println("\n---------17---------\n",V(rg(1,10)).fmap(e->e.ffx(p->p.real())));
        println("\n---------18---------\n",C(123).ifx(e->e*2).ifx(e->e*34));
        println("\n---------19---------\n",C(3,4).abs());
    }
}